/*
 *  Pacala pairwise.c
 *  
 *
 *  Created by Ben on 13/11/09.
 *  Copyright 2009 University of Sydney. All rights reserved.
 *
	Variant on Localiser.c
	Compiled function for running through R.
	Takes point data and grid size and calculates pairwise distances between points in current and neighbouring cells
	Only keeps distances below a threshold value
	Outputs the reults as a .csv

 */

# include <stdio.h>
# include <math.h>
# define pi 3.141593

void localiser(double *iX, double *iY, double *oX, double *oY, double *cellsize, int *popsize, int *oppsize, int *spaceX, int *spaceY, int *indexes, int *indID, double *dists, int *L);
void neighbours(int gX, int gY, int maxX, int maxY);


//Global variables
int n[9][2];

//Function to calculate pairwise distances between inds and opps.  
//Number of pairwise distance calculations is reduced by using Pacala grid
void localiser(double *iX, double *iY, double *oX, double *oY, double *cellsize, int *popsize, int *oppsize, int *spaceX, int *spaceY, int *indexes, int *indID, double *dists, int *L){
	
	double cs=*cellsize;
	const int gridX=(*spaceX)/cs, gridY=(*spaceY)/cs, maxopps=*oppsize/(gridX*gridY);
	int ocells[gridX][gridY], ogridX, ogridY, i, k=0;
	int test[gridX*gridY][3*maxopps];
	
	//zero cells matrix
	for (i=0; i<gridX; i++) {
		for (int j=0; j<gridY; j++){
			ocells[i][j]=0;
		}
	}
	
	// Assign opportunities (by their array index) into the test matrix and count opps per cell
	for (i=0; i<*oppsize; i++){
		ogridX=oX[i]/cs;
		ogridY=oY[i]/cs;
		test[gridX*ogridY+ogridX][ocells[ogridX][ogridY]]=i;
		ocells[ogridX][ogridY]++;
	}
	
	
	//Step through individuals and grab opps in each ind's cell plus neighbours 
	for (i=0; i<*popsize; i++){
		ogridX=iX[i]/cs;
		ogridY=iY[i]/cs;
		neighbours(ogridX, ogridY, gridX-1, gridY-1);
		int m=0;
		while (n[m][0]!=-1 && m<9) {
			int j=0;
			while (j<ocells[n[m][0]][n[m][1]]) {
				indexes[k]=test[gridX*n[m][1]+n[m][0]][j];
				indID[k]=i;
				j++;
				k++;
			}
			m++;
		}
	}
	//at this point, we should have an array of opp indexes and an array of ind indexes.  k=length of data in these arrays
	*L=k;
	
	//Calculate pairwise distances 
	for (i=0; i<k; i++) {
		dists[i]=pow(pow((iX[indID[i]]-oX[indexes[i]]), 2)+pow((iY[indID[i]]-oY[indexes[i]]), 2), 0.5);
		//temp[i]=kernel(dists[i], *u, *v);
	}
	
	//output results
	FILE *output;
	output=fopen("pdists.csv", "w");
	printf("Source , Potential , Distance \n");
	char a, b;
	for (i=0; i<k; i++){
		a=indID[i];
		b=indexes[i];
		printf(a, " , ", b, " , ", dists[i], "\n");
	}
	fclose(output);
}


//Function to return a matrix of cbind(X, Y) neighbouring cells.  Bounded space!
void neighbours(int gX, int gY, int maxX, int maxY){//gx, gy are grid cell refs (counting from 0), maxX is max value of gX
	
	int i, j, k, m, tempY, tempX;
	const int xpos[9]={-1, 0, 1, -1, 0, 1, -1, 0, 1}, ypos[9]={1, 1, 1, 0, 0, 0, -1, -1, -1};
	
	//Clear neighbours matrix
	for (i=0; i<9; i++){
		for (j=0; j<2; j++){
			n[i][j]=-1;
		}
	}
	
	// Walk down through 9 grid cells
	m=0;
	for (k=0; k<9; k++) {
		tempX=gX+xpos[k];
		tempY=gY+ypos[k];
		if(tempX>=0 && tempX<=maxX && tempY>=0 && tempY<=maxY){
			n[m][0]=tempX;
			n[m][1]=tempY;
			m++;
		}
	}
}

/*double kernel(double dist, double U, double V){
	//Fingers crossed I got this right!
	return 2;//(pow(u, v)*v*pow(pow(v, v)*pow(pow(u, 2)*v+pow(dist, 2), -2-v), 0.5))/(2*pi);
}*/


